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Abstract 

We construct asymptotic expansions for ordinary differential equa¬ 
tions with highly oscillatory forcing terms, focussing on the case of 
multiple, non-commensurate frequencies. We derive an asymptotic ex¬ 
pansion in inverse powers of the oscillatory parameter and use its trun¬ 
cation as an exceedingly effective means to discretize the differential 
equation in question. Numerical examples illustrate the effectiveness 
of the method. 
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1 Introduction 


Our concern in this paper is with the numerical solution of highly oscil¬ 
latory ordinary differential equations of the form 

M 

y'{t) = f{y{t)) + t > 0, y{0) = y^ e (1.1) 

m=l 

where f : ^ and ai, • • • , aM '■ K-i- —>• are analytic and ui,U 2 , -'' > G 

M \ {0} are given frequencies. We assume that at least some of these fre¬ 
quencies are large, thereby causing the solution to oscillate and rendering 
numerical discretization of (1.1) by classical methods expensive and ineffi¬ 
cient. Many phenomena in engineer and physics are described by the oscil- 
latoryly differential equations(Chedjou2001, Fodjouong2007, Slight2008 and 
so on). 

A special case of (1.1) with W 2 m-i = fnio, C02m = —moj, m = 0,1, ■ ■ ■ , [M/2J , 
where w S> 1, is a special case of 

OO 

y'{i) = fiyit)) + G{y) Y t > 0, y(0) = yo £ &, (1.2) 

k=—oo 

where G : x C'^ ^ C'^ is smooth, which has been already analysed at 

some length in (Condon, Deaho and Iserles 2010). It has been proved that 
the solution of (1.2) can be expanded asymptotically in 

OO ^ OO 

I/(()~Po.o(«) + E^ E <>0, (1.3) 

r=l m=—oo 

where the functions which are independent of uj, can be derived recur¬ 
sively; Pj. Q by solving a non-oscillatory ODE and ^ / 0, by recursion. 

An alternative approach, based upon the Heterogeneous Multiscale Method 
(E and Engquist 2003), is due to Sanz-Serna (2009). Although the theory 
in (Sanz-Serna 2009) is presented for a specific equation, it can be extended 
in a fairly transparent manner to (1.2) and, indeed, to (1.1). It produces 
the solution in the form 

OO 

y(t)~ (1.4) 

m=—oo 
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where Km{t) = 0(uj ^), m € Z. 0 Formally, (1.3) and (1.4) are linked by 


^m(^) — ^ 


E ^Pr,oi^), rn = 0, 

r=0 

cx> 

E l^Pr,mit), m / 0. 

I, r=l 


We adopt here the approach of (Condon et al. 2010), because it allows us 
to derive the expansion in a more explicit form. 

The highly oscillatory term in (1.2) is periodic in tui: the main difference 
with our model (1.1) is that we allow the more general setting of almost 
periodic terms (Besicovitch 1932). It is justified by important applications, 
not least in the modelling of nonlinear circuits (Giannini and Leuzzi 2004, 
Ramirez, Suarez, Lizarraga and Collantes 2010). 

Another difference is that we allow in (1.1) only a hnite number of dis¬ 
tinct frequencies in the forcing term. This is intended to prevent the oc¬ 
currence of small denominators, familiar from asymptotic theory (Verhulst 
1990). Note that (Chartier, Murua and Sanz-Serna 2012) (cf. also (Chartier, 
Murua and Sanz-Serna 2010)) employs similar formalism - a finite number 
of multiple, noncommensurate frequencies - except that it does so within 
the ‘body’ of the differential operator, rather than in the forcing term. 

We commence our analysis by letting Uq = {1,2, ■ ■ ■ , M} and ujj = KjUi, 
j = 1, • • • ,M, where a; is a large number which will serve as our asymptotic 
parameter. Consequently, we can rewrite (1.1) in a form that emphasises 
the similarities and identihes the differences with (1.2), 


y'{t) = f{y{t)) + ^ t > 0, y(0) = Vq ^ (1.5) 

m&lAa 


Section 2 is devoted to a ‘warm up exercise’, an asymptotic expansion of a 
linear version of (1.5), namely 


y'{t) = Ay{t) + ^ t > 0, y{0) = yo G (1.6) 

mgWo 


where A is a d x d matrix. Of course, the solution of (1.5) can be written 
explicitly, but this provides little insight into the real size of different com¬ 
ponents. The asymptotic expansion is considerably more illuminating, as 

^In the special case considered in (Sanz-Serna 2009) it is true that Km{t) = Ofa;”'"*'), 
m £ X, but this does not generalise to (1.2). 
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well as hinting at the general pattern which we might expect once we turn 
our gaze to the nonlinear equation (1.5). 

For the ODE with the highly oscillatory forcing terms with multiple 
frequencies, the asymptotic method is superior to the standard numerical 
methods. With less computational expense, this asymptotic method can 
obtain the higher accuracy. Especially, the asymptotic expansion with a 
fixed number of terms becomes more accurate when increasing the oscillator 
parameter w. 

In Section 3 we will demonstrate the existence of sets 


and of a mapping 

CX> 

a : Z/Zj. —)• M 

r =0 

such that the solution of (1.5) can be written in the form 

y{t) ~ Po,o(i) + t > 0. (1.7) 

‘ ^ UJ ' ^ 

r=l m£lAr 

As can be expected, the original parameters {ki, K 2 , • • • , km} form a subset 
of Ur- However, we will see in the sequel that the set a (Ur) is substantially 
larger for r > 3. In the sequel we refer to elements of aiJAr) as {am '■ fn G Ur}. 

The functions which are all independent of w, are constructed 

explicitly in a recursive manner. We will demonstrate that the sets Ur are 
composed of n-tuples of nonnegative integers. 

In Section 4 we accompany our narrative by a number of computational 
results. Setting the error functions are represented as 

es{t,uj) = y{t) - pQ Q{t) Pr■,m(^)e*'""*“^ 

‘ ^ OJ ' ^ 

r=l m£Ur 

we plot the error functions in the figures to illustrate the theoretical analysis. 
The expansion solvers are convergent asymptotically. That is, for every 
e > 0, fixed s, the bounded interval for t, there exists wq > 0 such that for 
CO > coq, the error function |es(t,a;)| < s. However, for increasing s, fixed t 
and 00 , we will come to the convergence of the expansion in future paper. 
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2 The linear case 


Our concern in this section is with the linear highly oscillatory ODE 
(1.6), which we recall for convenience, 


y' = Ay + ^ t > 0, y(0) = yo e (2.1) 


m&Uo 


Its closed-form solution can be derived at once from standard variation of 
constants. 



( 2 . 2 ) 


However, the finer structure of the solution is not apparent from (2.2) with¬ 
out some extra work. Each of the integrals hides an entire hierarchy of 
scales, and this becomes apparent once we expand them asymptotically. 

The asymptotic expansion of integrals with simple exponential operators 
is well known: given g G t) and \g\ S> 1, 



(Iserles, Nprsett and Olver 2006). Eor any m G Uq we thus take g{x) = 


e and y = (recall that Km / 0), therefore 

y{t) ~ 






We deduce the expansion (1.7), with him = {0}U7/o = {0,1, • • • , M}, m G N, 
and the coefficients 


Po,o(i) = 


5 













We thus recover an expansion of the form (1.7), where (Tq = 0 and 
o'm = fn G Uq. Note that linearity ‘locks’ frequencies: each Pj-m depends 
just on ttm, m ^Uq. This is no longer true in a nonlinear setting. 

3 The asymptotic expansion 

3.1 The recurrence relations 

We are concerned with expanding asymptotically the solution of 


y' = f{y) + t > 0, y(0) = y^ e (3.1) 


mdUo 


The function / is analytic, thus for every n G N there exists the nth differen- 


n times 


/ s 

tial of /, a function /„ : x C'^ x —)• C'^ snch that for every sufficiently 

small |t| > 0 



n=l 


Note that is linear in all its arguments in the square brackets. 
We substitute (1.7) in both sides of (3.1) and expand about 



m£Ui 




n=l 


£l=l ki£Ui^ 
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Pf k ^ 




A—/ ^in 

^n —1 

CO ^ OD OO 

/(po.o) + X 3 -i 'E^'"Y^ 


+ E 


CLm^ 


IKrn^t 


m^Uo 

1 


fl\ z-/ z-/ ^;^1+^2H-1-^' 

Ti=\ ^-^ = \ ^j^ = \ 


E ■■■ E ^'n(Po.o) [PilA- 


.p<.,*.]®‘‘"*-+-‘""*"’“‘+ E “'»" 


XKyTji^Ujt 


m^Uo 


cxD 1 cx) 1 


f^Po,o) + Y.-^Y.^ Y. Y --- Y /n(Po,o) [P£i,fcir-- ,P£„ 
n=l ■ r=n £el°,r fclSW^^ k„&Ui„ 


X 


giK^+-+a;,„M ^ 


ihvffiLiJt- 


m£Uo 


1 1 _^ ^ ^ 

fiPo,o) + Yz;^Y:^_ Y Y--- Y fn{Po,o)[PhM^--- ^Pe. 

r=i n=i ' £el°^^k-LeUe^ k„eUe^ 


xe*Ki+-+^'=nM+ amc 
m^Uo 


IK/ffiLOt’ 


where 

= {£ G N" : = r}, 1 < n < r. 

There is a measure of redundancy in the last expression: for example 
there are two terms in I 2 31 namely ( 1 , 2 ) and ( 2 , 1 ), but they produce identi¬ 
cal expressions. Consequently, we may lump them together, paying careful 
attention to their multiplicity. More formally, we let 


In,r = {■£ £ N"" : = r,ii < £2 < ■ ■ ■ < in}, 1 < n < r, 


the set of ordered partitions of r into n natural numbers and allow d£ stand 
for the multiplicity of £, i.e. the number of terms in 1 ° ^ that can be brought 
to it by permutations. For example, 0i^2 = 2, while for r = 4 there are five 
terms. 


04 — 1, 01,3 — 2, 02,2 — 1, 01,1,2 — 3, 01,1,1,1 — 1- 

We introduce the multiplicity and obtain a more compact form for the equa¬ 
tion 

Po,o + Y (3.2) 

mgWi 




+E 


r=l 


W 


E P'r,r 




m£Ur 
oo 


mdUr+l 


'L0'yY} (jjf^ 


f(Po,o) + Y.Zf'i2^_ Y. Y ••• /n(Po,o) ,P£„ 

2 ^ — 1 TL=\ ■^^^n,r kn^lA^^ 


X 


m^Uo 


IKi'ffiUjt. 


We separate different powers of uj in (3.2). The outcome is 

Po,0 + Y = f{Po,o) + Y 

mdlAi mdUo 

for r = 0 and 


(3.3) 


^ iamPr+i,me^‘'^"" (3.4) 

m^lAr m£Ur+i 

= Y^_Y^^ Y ••• /n(P 0 , 0 ) [P 4 ,fcl’--- >P^n,fcn] 

n=l ^Eln.r 

for r G N. 

Before we embark on detailed examination of the cases r = 0,1,2, fol¬ 
lowed by the general case, we must impose an additional set of conditions 
on the coefficients p^m- Similarly to the expansion of (1.2) in (Condon et 
al. 2010), we obtain non-oscillatory differential equations for the coefficients 
Pj.Q, r G Z+, which require initial conditions. We do so by imposing the 
original initial condition from (3.1) on pg q and requiring that the terms at 
origin sum to zero at every oj scale. In other words, 

Po,o(0)=l/0> Pr,o(0) = - Y1 Pr-,m(0), ^ G N. (3.5) 

m£Wr\{0} 


3.2 The first few values of r 

The expansion (1.7) exhibits two distinct hierarchies of scales: both 
amplitudes uj~^ for r G and, for each r G N, frequencies In 

(3.3) and (3.4) we have already separated amplitudes. Next we separate 
frequencies. 
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For r = 0 (3.3) and (3.5) yield the original ODE (3.1) without a forcing 
term, 

Po,o = fiPofi), t > 0, Po,o(0) = yo> 

as well as the recursions 

Pl,m = - -: m= I,-- - ,M 

IKm 

(recall that Km / 0). Therefore am = Hm, fn = I,-- - ,M. We set, for 
reasons that will become apparent in the sequel, 

= {0}UZYo = {0,l,--- ,M}, 


with (To = 0. 

For r = 1 we have lip = {1}, 9i = 1, and (3.4) yields 

E + E = E /i(Po,o) [Pi J 

mdUi mgW2 m&Ai 

We set 

U2=Ui = {Q,l,--- ,M} 

and (recalling that Pi m 3-re already known for m ^ 0) 

P'i,o = /i(Po,o)[Pi,o]> t > 0, Pi,o(0) = - Pl,m(0), 

meWi\{0} 

P2,m = 7^{/l(P0,0)[Pl,m] “ ^ G \ {0}. 

tKjn 

An explanation is in order with regard to our imposition of 0 G Z^i. We 
could have accounted for all the r = 1 terms in (3.4) without any need 
of the Pi 0 term. However, in that case the outcome would not have been 
consistent with the initial condition (3.5) and this is the rationale for the 
addition of this term. 

Our next case is r = 2. The case is not so straightforward to deduce. 
Since lip = {2} and l 2 p = {(1,1)}, we have from (3.4) 

E PUe“"""‘+ E = E 

mgW2 mgWa mgW2 

+ ^ E E /2(Po,o)bi.™,.Pi,„.Je‘''-‘+‘’""’“‘. (3.6) 

miSWi m2&Ai 
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We need to choose to match frequencies in the above formula. The set 
U 2 accounts for the frequencies 0 , ki, K 2 , •'' > but we must also account 
for Kj + Kj for i, j = 1, 2, • • • , M. Therefore we let 

^3 = ^2 U {(mi, m 2 ) : 0 < mi < m 2 < M}. 

We note two important points. Firstly, for i 7 ^ j, Ki + Kj can be obtained 
for (j, i), as well as for (i,j). Secondly, it might well happen that there 
exist i, j, /c € {1, • • • , M}, i < j, such that Ki + Kj = Kk -in that case we do 
not include {i,j) in This motivates the definition of the multiplicity of 
m ^ hl^\hl 2 (which we will generalise in the sequel to all sets Ur). Thus, 
for every < ii < £2 ^ M we let equal the number of cases when 

^ 7 r(£i) + ^ 7 r(£ 2 ) ~ where 7 r(€) is a permutation of i. Likewise, we let 
where 0 < ^1 < .^2 < Tf and 1 < mi < m 2 < M, be the number of 
permutations such that -|- ^^(£ 2 ) = Kmi + «^m 2 - 

We can now separate frequencies. Firstly, the non-oscillatory term, cor¬ 
responding to (To = 0. It yields the non-oscillatory ODE 

P 2 ,o =/i(Po,o)IP 2 ,o] + 2 Ph,e2f 2{Po,o)lPi,ei^Pi,e2\^ 

+^.^ 9=0 

ii<e2 

whose initial condition, according to (3.5), is 

P 2 ,o( 0 ) = - P 2 .m( 0 ). 

meU2\{0} 

Secondly, we match all the terms in ^^ 2 \{ 0 }, and this results in the recurrence 

'il<-mP3,m = fliPofi) [P2,m] “ P2,m + ^ A^2(Po,o) [Pl,£i > ^ 1 ,^ 2 ]' 

h<i2 

Finally, we match the terms in U^ \U 2 . Recall that these are pairs [mi, m 2 ) 
such that mi < m 2 and Km^ + Km^ / (Jj for j = 0,1, • • • , M. We obtain the 
recurrence 

i(Kmi +Km2)R3,(mi,m2) = ^ Y ^2 (Po,o) [Pl,£w Pi' 

~^^rn2 

h<i2 

(There is no danger of dividing by zero since we have ensured that Kmi + 
i^m 2 / = 0.) This accounts for all the terms in (3.6). 
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3.3 The general case r > 1 

We consider the ‘level r’ equations (3.4) noting that, by induction, the 
sets lAg are known for .£ = 1, 2, • • • , r and 0 € Ur- Moreover, we have already 
constructed all the functions ^s for m ^Ug, \ {0}, i = 1,2, - ■ ■ ,r, and all 
the PiqS for i = 0,1, • • • ,r — 1. The current task is to construct the set 
Ur+i, the functions Pr+i,m ^ Ur+i \ {0} and the function 

It will follow soon that all the terms in Ur+i are of the form Kj-^ + + 

• • • + where q < r and < j 2 < • • • < jq, We commence by setting 
number of distinct p-tuples (^1,^2; where 

- ■ ■ , "ig e {0,1, • • • ,M}, mi < m 2 <■■■ < mq, 

such that 

p q 

i=l i=l 

Examining the formula (3.4) we observe that the terms on the right hand 
side have exponents of the form where 

V = + (^k2 + • • • + O'ki £ ^£ii ) * = Ij ■ ■ ■ ) ^5 ^ £ Ir,n 

for some n € {1, 2, • • • , r}. It follows at once by induction on r that there 

exist g G {1, 2, • • • , r} and 0 < mi < m 2 < • • • < mq < M such that 

<? 

ri = 

i=l 

It might well be that such rj can be already accounted by Ur , in other words 
that there exists m £ Ur such that r] = am- Otherwise we add to Ur the 
ordered g-tuple {mi, m 2 -, ■ ■ ■ ,'mq). This process, applied to all the terms on 
the right of (3.4), produces the index set Ur+i 

Ur+i =Ur\^ {(w-i) • • • ) 0 < mi < m 2 < • • • < uig < M, g G 1, 2, • • • , r} . 

We impose natural partial ordering on Ur- first the singletons in lexico¬ 
graphic ordering, then the pairs in lexicographic ordering, then the triplets 
etc. This defines a relation mi < m 2 for all mi, m 2 G Ur- We let 

{ n q 

{(-, fc) . ki G U(^^, (- G IIr!,,r) y ^ ££ki — y ^ ^rrii ; kn, mi ^ ^ mq 

i=l i=l 
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where m GUr and n € {1, 2, • • • , r}. 

Let us commence our construction of recurrence relations by considering 
m G Ur \{0}. In that case we have 


io'mPr+l,m ~ Pr,m 
^ 1 

+ Y. Pfc/n(P0,0) ■ (3-7) 

n=i ■ t&n.T {e,k)ew^„^ 

Next we consider m G Ur+i \ Ur- Now the hrst sum on the left of (3.4) 
disappears and the outcome is 


^^mPr+l,m 


Y. ^fc/n(Po,o) ,P£„,fcJ 

n=i £ein,r 


(3.8) 


Finally we cater for the case m = 0: now the recurrence is a non-oscillatory 
ODE, 


Pr,0 


Y PkfniPofi)[P£i,ki^'"^Pir^,kn]^ 

n=i ■ £en„,r (£,fc)ew;^o 


t > 0, 
(3.9) 


Pr,o(9) ^ ^ Pr,m(9)' 

m£Ur\{0} 

For example, in the case r = 3 we have 
Il ,3 = { 3 }, l2,3 = {( 1 , 2 )}, 13^3 = {( 1 , 1 , 1 )}, 6*3 = ^*1,1,1 = Ij ^*1,2 = 2 , 

while 

U3 = { 0 , 1 , • • • ,M}U{(mi,m2) : mi < m2,Kmi+Km2 / Vm = 0 , • • • ,M}. 
We thus deduce from ( 3 . 7 ) that 

il^mPA,m = -p'3,m + /l (Po,o) [P3.m] + Y /2 (Po.o) [Pl,ji > ^ 2 , 12 ] 

jl<j2 

+ l Y PhJ2,nf3iPo,o)[Pi,n^Pi,j2^Pij3] 

^31 '^^32 m 

jl<j2<j3 
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for all m G ZY 2 \ {0} and 

i{^mi “t“ ^rn,2 )p4,(mi,7722) ^3,(mi,m2) (Po,o) [i^3,(mi,m2)] 

+ pTijrf2iPo,o)[Pi,n.P2,j2] 


—^mi G^m2 

h<h 


6 


'tjl G^j 2 —^TTii G^m2 

jl<j2<j3 


for {mi, m 2 ) G \ U 2 - Next we use (3.8): 

+ ^^^2 + ^"33)^4,(mi,m2,m3) 

1 
6 


/s iPo,o)[Pi,n ’ Pi,i2 > Pija] 


ii<j2<i3 

for all 1 < mi < m 2 < m 3 < M such that Kmi + «^m 2 + ^m 3 7 ^ c^m for 
m 

Finally, we invoke (3.9) to derive a non-oscillatory ODE for q, namely 
P3,0 = /l(Po,o)[P3,o] + Y1 Ph,j2f2{Po,o)\Pl,n,Pl,j2] 

^3\ ~^^32 

+ g Z5jl,i2j3'^3(PO,o)[Pl,Jl)Pl,j2’Pl,J3]’ 

jl<j2<j3 

P 3 ,o( 0 ) = - P 3 ,m( 0 )- 

meWaVto} 

3.4 A worked-out example 

Let M = 3 and 

KI = 1, K2 = V 2 , K3 = —1 — \/2. 

Therefore Ui = IA 2 = {0,1,2, 3}, 

Uq = 0, CTi = 1, (72 = \/ 2 , (73 = -1 - \/2 

and 

Pk — ^k,mi k,m — 0,1, 2,3. 
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Consequently, = /(po,o)> Po,o( 0 ) = 2 /( 0 ) and 


ai 02 

Pl,3- 


as 


(1 +V2)i' 

We commence with r = 1 . The ODE is now 

p'lfi = /i(Po,o)bi,o]> t > 0, Pi,o(0) = -Pi,i(0) - Pi,2(0) - Pi,3(0), 
while the recurrences are 

P2,l = -{/l(Po,o)[Pl,l] -Pl,l}’ P2,2 = ^{/l(Po,o)[Pl,2] -Pl,2} 

1 V zz 


1 


-{/l(Po,o)[Pl, 3 ] -p'l.s}- 


(1 + ^)i 

Note that the first two ’’levels” of the expansion are 

»(*) » Po,o + - [pi.o + + Pi,2e‘'^"‘ + Pi,3e-‘l‘+'^>"‘ 

Uu L 

Next ioU^ = {0,1,2, 3, (1,1), (1, 2), (1,3), (2, 2), (2, 3), (3, 3)}, with oi,i = 
2, (Ti ,2 = 1 + \/2, ci' 1,3 = —\/2, (T 2,2 = 2\/2, ( 72,3 = —1 and 173,3 = —2 — 2\/2. 

The only way to obtain uo = 0 using two terms from lAi is 0 + 0, therefore 
Pq q = 1, otherwise = 0. However, to obtain am for m = 1, 2,3 we have 
two options: 0 + m and m + 0. Therefore p^^ = 2, otherwise P^ = 0. For 

mi,m 2 i. J.U p Id 2,2 3,3 -I 1,2 1,3 2,3 r) ,i 

PeiX Pi,I = Pi 2 = P 3,3 = Pi ,2 = Pi, 3 = P 2,3 = 2, otherwise 

the coefficient is zero. Therefore 

P2,0 = /l(Po,o)[P2,o] + 2'^2 (Po,o)[Pi,05Pi,o]) 

P2,o(0) = -P2,l(0) -P2,2(0) -P2,3(0)> 

and the 0{uj~‘^) terms are 

^b2,o +P2,ie'"‘ +P2,2e‘'^‘ +P2,3e-‘I‘+'^>"‘1. 


Moreover, 


P 3 ,l = -{/i(Po,o)IP2,i] -P2,l + /2(Po,o)[Pl,0>Pl,l]}> 
P 3,2 = ;^{/l(Po,o)[P2,2] - P'2,2 + /2(Po,o)[Pi,0>Pi,2]}2 
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1 


P3,3 = 
P3,(l,l) 

P3,(l,2) 

P3,(l,3) 

P3,(2,2) 

P3,(2,3) 

P3,(3,3) 


:{/l(Po,o)[P2,3] -P2,3 + /2(P0,0)[Pl,0>Pl,3]}> 


(1 + V2)i 

= -^■f2(Po,o)[Pl,l,Pl,l], 


~ 7 r~i 7 ^/ 2 (Po,o)[Pl,15^1,2]) 

(1 + v 2 )^ 

= 2iPo,o)[Pl,l^Pl,3]i 

~ 4 v^^'^ 2 (Po,o)[Pi, 2 )Pi, 2 ]! 

= “■^/2(Po,o)[Pi,2)Pi,3]5 


The number of terms in is significantly larger: we need to add to Us 
a further nine terms - altogether we have 19 terms, which are displayed in 
Table 1. All the nonzero are displayed there as well. Note that there 

is no (1,2,3) term, because ki + ^2 + ^3 = 0, hence it is counted together 
with zero. 

In particular, 

P3,0 = /(Po,o)[P3,o] + /2(Po,o)[Pl, 01 ^ 2 , 0 ] + 3iPo,o)lPl,O^Pl,O^Plfi]-! 

3 3 3 

P 3 ,o( 0 ) = - 5 ^P 3 j( 0 ) - ^I^P 3 ,(i,£)( 0 )- 
i=i i=i i=j 


We conclude that the 0(u! terms in the asymptotic expansion of y are 


1 


Aujt 


P 3,0 + P 3 ,l^''~“ + P 3,2 

i{l-\-y/2)u)t 


U 

+P3,( 1,2)6 


+ P3 + P3,(1,1)6 

n 0-\ fOi. i-f 

P3,(2,2)6 


i2ujt 


n „-iV2uit I i2y/2u)t I 

+ P3,(1,3)6 +P3,(2,2)6 +P3,(2,3)6 


I „ ^-i{2+2V2)u)t 

+P3,(3,3)6 

All this can, of course, be carried forward to larger values of r. 


3.5 Two non-commensurate frequencies 

An interesting special case is M = 2, where, without loss of generality, 
7 ^ 0 is rational and K 2 irrational: this means that the only integer solution 
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Table 1: Ordered elements of 


m 


pT 

0 

0 

II 

0 

/^0,0,0 = l,Pl,2,3 = 2,P2,1,3 = 2,P3p 2 = 2, 

1 

CJl = 1 

/^o,o,i = 3, 

2 

(72 = \/2 

/^0,0,2 = 3, 

3 

1 

I 

II 

CO 

b 

Po,0,3 — 3, 

(1,1) 

(Tip — 2 

1,1 o 

Po,l,l ~ 

(1,2) 

(Tip = 1 + ^/2 

1,2 e 

P0,l,2 — 6, 

(1,3) 

O'lp = — 

1,3 

/^0,1,3 ~ 

(2,2) 

(T2p = 2\/2 

2,2 0 

P0,2,2 — "J, 

(2,3) 

0'2,3 = —1 

2,3 

P0,2,3 — '3, 

(3,3) 

(T3p = —2 — 2\/2 

3,3 0 

^0,3,3 ~ "J, 

(1,1,1) 

O'!, 1,1 = 3 

1,1,1 1 

Pl,l,l ~ 1, 

(1,1,2) 

o'lpp = 2 + \/2 

1,1,2 q 

^•1,1,2 ~ "J, 

(1,1,3) 

o'lpp = 1 — 

1,1,3 q 

Pi,1,3 ~ 

(1,2,2) 

0'1,2,2 = 1 + 2-\/2 

1,2,2 q 

Pl,2,2 ~ 

(1,3,3) 

o'lpp = —1 — 2^/2 

1,3,3 q 

Pl,3,3 ~ "J, 

(2,2,2) 

<^2,2,2 = 3\/2 

2,2,2 , 

P2,2,2 ~ ^1 

(2,2,3) 

0'2,2,3 = —1 + 

2,2,3 q 

P2,2,3 ~ "J’ 

(2,3,3) 

0'2,3,3 = —2 — \/2 

2,3,3 q 

P2,3,3 ~ '^5 

(3,3,3) 

0'3,3,3 = —3 — 3-\/2 

3,3,3 

P3,3,3 ~ 


to miKi + m 2 K 2 = 0 is ml = m2 = 0. This simplifies the argument a great 
deal. 

Simple calculation now confirms that 
P'ofl = /(Po,o)> Po,o(0) = Po> 

Pl,Tn 5 ^ 

I Km 

p'i,o = /i(Po,o)[Pi,o]> t > 0, Pi,o(0) = -Pi,i(0) - Pp2(0), 

P2,m = 7^{/l(P0,0)IPl,m] “ m = 1, 2, 

Ir^m 

P' 2,0 = /l(Po.o)[P 2 ,o] + ^/ 2 (P 0 , 0 )[Pl, 0 >Pl, 0 ]> P2,o(0) = -P2,l(0) “ P2,2(0)> 

P3,m = T^{/l(P0,0)[P2,m] -P2,m + /2(Po,o) [Pl,0> Pl,m]}> m = 1, 2, 
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1 


P3,{m,m) f 2(.Po,o)[Pl,mi Pl,m\'! ^ 1) 

P3,(1,2) = -^;^^^f2iP0,0)[Pl,l,Pl,2], 


P3,0 = /l(Po,o)[P3,o] + f 2 iPQ,o)\Pl,OjP 2 ,o\ + g /3(Po,o) [Pl, 0 ) Pl, 0 ) Pl,o]) 

P3,o(0) = “P3,l(0) “ ^ 3 , 2 ( 0 ) “ P3,(l,l)(0) “ P3,(l,2)(0) “ P3,(2,2)(0)' 
The next ‘generation’ is 


= — {/l(P0,0)[P3,m] -P3,m + /2 (Po,o) [Pl,0> P2,m] 

IKm 

+ f2iPo,o)[Pl,m,P2,o] + ^/3(P0,0)[Pl,0>Pl,0>Pl,0]| > m = 1,2, 
P4:,{m,m) 2>iK "1*^ 1 (^^0,o) [i^3,(m,m)] f 2i.Po,o'}[Pl,rn^ P2^m\ 

+ ^/3(P0,0)[Pl,0>Pl,m>Pl,m]| > m = 1,2, 

P4,{1,2) = _|_ ^2) {•^l(Po>o)[P3,(l,2)] + /2(Po,o)[Pl,l)P2,2] 

+ /2(Po,o)[Pl,2)P2,l] + /3(Po,o)[Pl,0)Pl,l)Pl,2]} ) 

P4,{m,m,m) ^ ^ iPo,o)\Pl^m^ Pl,m'> Pi,ml 5 ^ 

^^’(^’2,2) = ^^^^^^^^/3(Po,o)[Pi,1>Pi,2,Pi,2] 


and 

P4,0 = /l(Po,o)[P4,o] + /2 (Po,o)[Pi,05P3,o] + 2-^2 (Po,o) [^2,05 P 2 ,o] 

+ ^/3(P0,0)bl,0>Pl,0>P2,0] + ^/4(P0,0)bl,0>Pl,0>Pl,0>Pl,0]> 

P4,o(0) = “P4,l(0) “^ 4 , 2 ( 0 ) “P4,(l,l)(0) “P4,(l,2)(0) “P4,(2,2)(0) 
“P4,(l,l,l)(*^) “ P 4 ,(1,1,2) (0) “ P4,(l,2,2)(0) “ P4,(2,2,2)(0)- 

Greater, bnt not insnr mount able effort is required to develop a general 
asymptotic expansion in this case. However, to all intents and purposes, 
expanding up to r = 4 is sufficient to derive an exceedingly accurate solu¬ 
tion. 
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3.6 Comments 


Comment 1: As we increase levels r, we are increasingly likely to en¬ 
counter the well-known phenomenon of small denominators (Verhulst 1990): 
unless Krn = m = 1, • • • ,M, where all the 'ipmS are rational (in which 
case, replacing w by its product with the least common denominator of the 
'ipmS, we are back to the case (1.2) of frequencies being integer multiples of 
ui) and positive, the set of all finite-length linear combinations of the k^s is 
dense in M (Besicovitch 1932). In particular, we can approach 0 arbitrarily 
close by such linear combinations. This is different from KmS summing up 
exactly to zero: as demonstrated in Subsection 3.4, we can deal with the 
latter problem but not with the denominators in (3.7) or (3.8) becoming 
arbitrarily small in magnitude. Like with other averaging techniques, there 
is no simple remedy to this phenomenon. This restricts the range of r at 
which the asymptotic expansion is effective. Having said so, and bearing 
in mind that the truncation of (1.7) to r < ii yields an error of 
and, |a;| being large, we are likely to obtain very high accuracy before small 
denominators kick in. Hence, this phenomenon has little practical implica¬ 
tions. 

Incidentally, this is precisely the reason for the requirement that, unlike 
in (1.2), the number of initial frequencies is finite. Otherwise, we could 
have encountered small denominators already for r = 3 and this would have 
definitely placed genuine restrictions on the applicability of our approach. 

Comment 2: There is an alternative to our expansion. We may decide 
that the are symbols, rather than specific numbers. Not being assigned 
specific values, it is meaningless to talk about sets because ki^ + 

K 12 = say, has no meaning (except when £i = 0, £2 = m). Of course, 
in that case we may have several distinct terms which correspond to the 
same frequency, once we allocate specific values to the KmS, but the quid 
pro quo is considerable simplification and no multiplicities (which depend 
on specific values of k^s, hence need be re-evaluated each time we have 
new frequencies). Unfortunately, this approach has another, more critical, 
shortcoming. We must identify all linear combinations of that sum up 
to zero, because they require an altogether different treatment. 

It would have been possible to proceed differently, by separating all the 
terms in of the form subsets: those that sum to zero and 
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those that are nonzero. We do not need to specify which is which - this 
becomes apparent only once values are allocated - just to remember that 
the nonzero sums give rise to new frequencies, with coefficients derived by 
recursion, while zero sums are lumped into a differential equation for a non- 
oscillatory term. While this is certainly feasible, it seems that the current 
approach is probably simpler and more transparent. 

4 Numerical experiments 

In the current section we present two examples that illustrate the con¬ 
struction of our expansions and demonstrate the effectiveness of our ap¬ 
proach. In each case we compare the pointwise error incurred by a truncated 
expansion (1.7) with either the exact solution or the Maple routine rkf45 
with exceedingly high error tolerance, using 20 significant decimal digits. 
Specifically, we measure the components of 

vit) - Po,o{*) - E E 

r=l mSWr- 

for different values of s. 

4.1 A linear example 

We consider the equation 

x + -x + —x = te^^‘^^+ t>0, x(0) = x(0) = 0.5. (4.1) 

5 5 

Letting y = [x,x]"^, we reformulate (4.1) as the system 

.f , t>0, y(0) = 

This being a linear equation, the exact solution and its asymptotic expansion 
are available explicitly using the theory from Section 2. 

Figures 4.1 and 4.2 display the real part of the error functions in com¬ 
puting X and X , respectively, for s between 0 and 4 within t G [0,5]. It is 
clear that each time we increase s, the error indeed decreases substantially, 
in line with our theory. 
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tOlh-HOl 







1. x 10 
8.x 10 
6.x 10 
4.x 10 

2. x 10 


-2.x 10 


-4.x 10 
-6.x 10 
-8.x 10 



Figure 4.1: The real part of the error in x committed by the asymptotic 
expansion, as applied to the linear system (4.1) with uj = 500 for s = 
0,1, 2,3,4 (from top left onwards). 
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Figure 4.2: The real part of the error in x committed by the asymptotic 
expansion, as applied to the linear system (4.1) with uj = 500 for s = 
0,1, 2,3,4 (from top left onwards). It is shown that the error decreases once 
s is increased from 0 to 1. 
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Identical information is reported in Figs 4.3 and 4.4 for frequency oj = 
5000 within t € [0,5]. A comparison with the two previous hgures empha¬ 
sises the important point that the efficiency of the asymptotic - numerical 
method grows with uj, while the cost is to all intents and purposes identi¬ 
cal. Indeed, wishing to produce similar error to our method with s = 4, 
the Maple routine rkf45 needs be applied with absolute and relative error 
tolerances of 10“^^ and 10“^® respectively. H Although the method is robust 
enough to produce correct magnitude of global errors, this comes at a steep 
price. Thus, while our method takes less than one second to compute the so¬ 
lution and requires ~ 17.5 kbytes of storage, rkf45 takes ~ 2740 seconds to 
compute the solution for oj = 500 and requires ~ 10^ kbytes. This increases 
to 4533 seconds and 1.6 x 10^ kbytes for uj = 5000. 


4.2 A nonlinear example in Memristor circuits 

In this subsection, the method in this paper is developed for Memristor 
circuits subject to high-frequency signals. Consider the following differential 
equation governing a circuit with two memristors similar to that given in 
BoCheng2011, 


= ysit) 

2 / 3 W = ay 3 it){d - (1 -k 3yf (t))) - 

_ iy3{t)-y4{t)){l+3y^{t)) 


a-b3[t)-y4,(t))il+3y^ (0) 

l+e(l-l-3y|(0) 


2/4 W = 


l+e{l+3yl (t)) 

y'S) = -byAit) - cy^it) + sit) 


+ 2/5(t) 

t > 0 


(4.2) 


with the initial conditions 


(yi(0),y2(0), 1 / 3 ( 0 ),y4(0),y5(0))^ = (ci,C2,0, 10 '^,0)'^ = yo 


where o = 8 , 6 = 10, c = 0, d = 2, e = 0.1, ci = —0.8, C 2 = —0.4. 
The unknown functions are yi (t ), 1/2 (t) j 2/3 (t) j 2/4 (t) > 2/5 (t) • The forcing term 
is sit) = ^ ^ ki = 1 , 

Ki = —1, K 3 = V 2 , K 4 = —\/2 and u is our oscillatory parameter. The circuit 

^The fact that Figs 4.1a and 4.1b are identical is a fluke-anyway, it is evident from the 
results on Page 5. 

^Such error tolerances are impossible in Matlab, which explains our use of Maple. 
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Figure 4.3: The real part of the error in x committed by the asymptotic 
expansion, as applied to the linear system (4.1) with u = 5000 for s = 
0,1, 2,3,4 (from top left onwards). 
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Figure 4.4: The real part of the error in x committed by the asymptotic 
expansion, as applied to the linear system (4.1) with u = 5000 for s = 
0,1, 2,3,4 (from top left onwards). 
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figure is shown in Fig. 4.5 where the corresponding parameter relationship 
between Fig. 4.5 and the equation (4.2) are 

4>i = 2/1, 4>2 = 2 / 2 , 113 = 2/3, 114 = 2/4, h = 2/5, ^2 = 1, a = 1/Ci, 

Wi = l + 2>yl, VF 2 = 1 + 31/2, b = l/L, c = rlL, d = G, e = R. 


- sit) + + ^2(0 - R 



Figure 4.5: The Memristor circuits. 


This circuit equation can be written in vector form 


y'(t) = f{y) + t > 0 , 

where y{t) = (i/i(t), 2/2(t), 2/3(i), 2/4(t), 2/5(i))^, 


f{y) 


( 2/3 

y4-y3 

l+e{l+Zyl) 

ayz{d- (1 + 3y?)) 

(y3-y4)(i+3y^) 
l+e(l+3y^) ^ 

\ -byi- cy 5 


\ 

a+3-y4)(l+3yi) 

i+e(T+%f) 




and 



f 

0 

\ 


f 

0 

\ 


f 

0 

\ 


f 

0 

\ 



0 




0 




0 




0 


ai{t) = 


0 


,a2it) = 


0 


,a3(t} = 


0 


,a4(t) = 


0 




0 




0 




0 




0 



V 

Ab 

2i 

/ 


V 

Ab 

2i 

/ 


V 

Ab 

2i 

/ 


V 

Ab 

2i 

/ 
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or in a more compact form as 

4 

y'{t) + f{y{t)) = Y, t > 0. 

m=l meUo 

(4.3) 

where Uq = {1,2, - ■ ■ , 4} is an initial set and ujj = kjuj, j = 1, 2, • • • . 
Then the asymptotic method is developed for this type of equation. 

The formation of the terms in the asymptotic method for the given 
memristor system shall now be described. 

4.2.1 The zeroth terms 

Denote {p)j as the j-th. element of the vector p. When r = 0, set 
Ui = {0,1, 2, 3,4}. Then the zeroth term Po o(^) obeys 

Po,o =/(Po,o)> Po,o(0) = I/O = (ci,C 2 ,0,10"^, 0)^, 


where 


/(Po,o) = 


/ (Po,o)3 \ 

(Po,o)4~(Po,o)3 

l+e(l+3(pg 0 ) 2 ) 

\ M ^2^^ “((Po,o)3-(Po,o)4)(l+3(Po,o)i) 

«(Po, 0 ) 3(4 - (1 + 3(Po,o)l))-i+e(l+3(po,o)^)- 

((PO,o)3~(Po,o)4)(1+3(Po_o)2) I ^ 

l+e(l+3(po,o)^) + 

V -^(Po,o)4 - C(P0,0)5 / 


In addition, the recursions enable the determination of Pi ^{t), m 7 ^ 0, 

\ 


Pi,lit) = — 

IKl 


/o \ 
0 
0 
0 

\ 2i / 


r ^ 1 

Pl,2it) = — 

IK 2 


( 0 
0 
0 
0 

v 2i / 





\ 

1 

0 

II 

1 

0 

0 

0 

m3 

0 

IK4 

0 


\§) 


v-f / 


The corresponding derivatives are p\ ^{t) = p'^ 2 (t) = p'l ^it) = p'l ^{t) = 0. 
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4.2.2 The r = 1 terms 

For r = 1, set IA 2 = U\ = {0,1, 2,3,4}. This yields 

P' 1,0 =/i(Po,o)bi,o]> 


/o \ 

0 
0 
0 

V f(2 + V2) j 

P2,m = 7^{/l(P0,0)[Pl,m] “ P'l,m} = (Po,o)[Pl,m]> m £ \ {0}, 


Pl,o( 0 ) = - Pl,m( 0 ) 

meWi\{0} 


IKn 


ZKr, 


where 


P2,i(i) = 77 


(zKl)" 


P2,3(^) = TT 




( 0 \ 
0 
0 

M 
2* - 

L —cAb . 
\ 2i / 

/o \ 
0 
0 

M 
2i , 

I —cAb ; 
\ 2i / 


P2,2{t) = -p 


{in2y 


P2,4(^) = T 




/O \ 
0 
0 

_M 

\ / 

\ 2i / 

( 0 \ 
0 
0 

2i 

{ <^^46 , 

\ 2i / 


4.2.3 when r = 2 

The r = 2 layer is the first layer in which additional frequencies must be 
considered and we set 


Us = {0,1,2,3,4, (1,1), (1,3), (1,4), (2,2), (2,3), (2,4), (3,3), (4,4)}. 

Note that the (1,2) term and (3,4) terms are not present as addition of these 
frequencies would result in zero which is present in the set. 

We will first consider the P 2 ,o(^) term. Since Pii = I, Pi 2 = ‘2. and 
P 3 4 = 2 , the term p 2 ,o satisfies 

P2,0 = /l(Po,o)[P2,o] + 2 P?i,£2'^2(Po,o)[Pi,£i)Pi,£2] 

K£, -\-K£^=0 
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— /l(Po,o)[P 2 ,o] + 2 '^ 2 (Po,o)[Pi, 0 )Pi,o] 

+ /2(Po,o)[Pl,l)Pl,2] + /2 (Po,o)[Pi,3!Pi,4] 

with the initial condition 


P2,o(0) = - Y1 P2,m(0) = 0, 
m£U 2 \{ 0 } 


where 


/2(P0,0)bl,m,Pl,fc] 


/ PlmMl{Pofi)Pl,k \ 
pTmM 2 {Po,o)Pl,k 
Pl,mM3{Po,o)Pl,k 
Pl,mM4{Po,o)Pl,k 
\ Pl,mM5iPo,o)Pl,k ) 


and Mj[y) is the 5x5 dimensional Jacobian matrix evaluated at the vector 
function Po,o(^) 


Mj{Po,o) 


/ d^fj d^fj d'^fj d^fj \ 


dyidyi 

d^fj 

dy2dyi 

d^i 

dyidy 2 

d^fi 

dyidyi 

d^fi 

dyidyi 

d^fj 

dyidya 

d^fi 

dyidyi 

d^fj 

dyidyi 

d^fi 

dyidys 

d^fj 

dyidys 

d^fi 

dysdyi 

dyidy2 

d^fi 

dyadys 

d^f, 

dyadyi 

d^fi 

dyadys 

d^fi 

dyidyi 

d^fi 

dy^dy2 
d fj 

dyidya 

d^fi 

dyidyi 
d fj 

dyidys 
d fj 

dysdyi 

dy5dy2 

dysdya 

dysdyi 

dysdys 


Furthermore, due to the fact that the first four elements of Pi^{t), 
m / 0, are zero, the nonoscillatory equation for P 2 fi{t) simplifies to 


P 2,0 = /l(Po,o)[P 2 ,o] + ^/ 2 (Po,o)[Pl, 0 >Pl,o]> P 2 ,o( 0 ) = 0, t > 0. 

Now consider the set U 3 and the recursions for P 3 First, we match all 
of the terms in U 2 \ {0} as these are also in U 3 , 


il^mP3,m = fliP0,0)[P2,m]-P2,m + l^ ^^,£2 /2 (Po,o) [Pl,£i > Pl,^] 


—^rn 

^1<^2 


= fl{P0,0)[P2,r 
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We then match to the remainder of the elements in lA-^. 






3,(mi ,m 2 ) 


\ /5“'-^2(Po.o)[Pmi>Pm2] = 0- 

+/t^2 -^m-i +^7712 

h<l 2 


Because of the nature of the elements of Pi ^ and P 2 ^ 
terms of P 3 ^(t), m / 0 , are 


the non-zero 


P3,i(^) = 


P3,2(i) = 


{iK.2f 


PsA*) = 


P-SAit) = 


( 0 


(iKl)3 


\ 


1 Ah 

l+e(l+3(po,o)2) 
a(l+3(Po,o)l) Ah 


2i 


l-|-e(l-|-3(pQ Q) 

— (l+SlPo.oli) Ab 

l+e(l-l-3(po,o)2) 
c^Ab 


+ 


—cAb 

2i 


Ab^ , 
2i 


2i 


( 0 


-Ab 


l+e(l+3(po_o)2) 
a(l-|-3(pQ o)i) —Ab 
l+e(l+3(Po,o)2) 

-(l+3(poo)|) cAb 

l-l-e(l-|-3(po,o)i) 

A62 c^Ab 
2i 


2i 


( 0 


(^ 3)3 


Ab 


1+®(i+3(Po,o)2) 
a(l+3(Po,o)l) Ab 
l+e(l+3(po,o)2) 

-(i+3(po^o) 2) Ab I -cAb 

l+e(l+3(Po,o)2) 


/ 


(iK4)3 


Ab^ I c^Ab 
2i 2i 

0 \ 

1 -Ab 

l-|-e(l-|-3(pg q)|) 2i 

a(l-|-3(pQ 0)2) —Ab 

l-|-e(l-|-3(po 0)2) 

-(l+3(Po,o)i) -Ab , cAb 
l-l-e(l-|-3(po,o)i) 

Ab^ _ e^Ab 
2i 2i 
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4.2.4 The r = 3 terms 


When r = 3, we note that / 9 qq = 1, p \2 = 2, = 2, /Oqo^o = 

Pop ,2 — 6 ) Po, 3,4 = 6 . Hence, the equation for P 3 q is 

P3,0 = /l(Po,o)[P3,o] X/ P\,i 2 f 2 {Po,o)[Pl,ii^P 2 ,l 2 ] 

=0 

^1<^2 

+ ^ Ph,t2,hhiPo,o)[Pi,ivPi,e2^Pi,h\ 

il<i2<^3 

= /l(Po,o)[P3,o] + /2(Po,o)[Pl,0)P2,o] + 2/2 (Po,o)[Pi,1)P2,2] 

+2/2(Po,o)[Pi, 3)P2,4] + g/3(Po,o)[Pl,05Pl,0)Pl,o] 

+ /3(P0,0)[Pl,0>Pl,nPl,2] + /3(P0,0)[Pl,0>Pl,3>Pl,4] 

= /l(P0,0)[P3,0] + /2(P0,0)[Pl,0>P2,0] + ^/3(P0,0)[Pl,0:Pl,0>Pl,0] 

+ /3(P0,0)[Pl,0>Pl,l>Pl,2] + /3(P0,0)[Pl,0>Pl,3>Pl,41- 


wit h 


P3,o(0) 


Y P3,m(0) 

meWsVfO} 


I 0 


Ab 

' l+e(l+3c|) 

A&a(l+3c2) 




l+e(l+3c|) 

cAb + ji^il±2£2) 1 (1 + ^ 

+ l+e(l+3c^)J V + 2V2 

(^ + 


\ 


/ 


where 


(/3(Po,o))i[Pi ,mi) Pi,1712 ’ Pi , 1713 ] 

5 5 5 d^f 

= YYY dyk,dyldyk. 

fcl = lfc2 = lfc3 = l ^ ^ ® 


)fcl(Pl )fc2(Pl , 1713 ) ks- 


Therefore, the asymptotic expansion including terms up to r = 3 is 

Pit) ~ Po,o{t) + - [Pi,o(t) +Pi, 2 (i)e*"^"* 

UJ 
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[P2,o(0 +P2,4We“^ 

+ ^ [P3,o(0 +P3,2(i)e‘"^"‘ +P3.3We*"^“* +P3,4We“^ • 

4.2.5 Numerical experiments 

The nonlinear Memristor circuits do not have a known analytical so¬ 
lution, we come to a reference solution, the Maple routine rkf45 with the 
accuracy tolerance AbsErr = 10“^^ and RelErr = 10“^^. The terms Pi q, 
P 2 0 and P 3 0 satisfy the non-oscillatory ODEs which is solved by the Maple 
routine rkf45 with AbsErr = 10“^® and RelErr = 10“^*^. 

Figures 4.6 to 4.10 show the error functions for yi, y 2 , 2 / 3 , 2/4 and 2/5 for 
the truncated parameter s = 0,1,2,3 within t G [0,3] when the oscillatory 
parameter is cu = 100 and oj = 1000. The error is seen to greatly reduce 
with an increasing number of r levels. Furthermore, with increasing the 
oscillatory parameter, the error of the asymptotic method decreases rapidly, 
a very important virtue of the method. 

In addition, the CPU time is compared to the Runge-Kutta method 
(rkf45) whose tolerance equals to 10“^^. It takes 275 seconds for uj = 500 
and about 2992 seconds for oj = 5000, respectively. The CPU time for the 
asymptotic method is about 11 seconds for oj = 500 and 13 seconds for oj = 
5000. As evident from our previous theoretical analysis, the computational 
cost is about the same for the asymptotic method regardless of the value of 
the oscillatory parameter. 
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Figure 4.6: The top row: the real parts of error function with s = 0 (the 
left) and s = 1 (the right) for yi with oj = 100. The middle row: the real 
parts of error function with s = 2 (the left) and s = 3 (the right) for yi with 
uj = 100. The third row: the real parts of error function with s = 0 (the 
left) and s = 1 (the right) for yi with uj = 1000. The fourth row: the real 
parts of error function with s = 2 (the left) and s = 3 (the right)for y-i with 


























































































































































Figure 4.7: The top row: the real parts of error function with s = 0 (the 

left) and s = 1 (the right) for y 2 with oj = 100. The middle row: the real 

parts of error function with s = 2 (the left) and s = 3 (the right) for y 2 with 
uj = 100. The third row: the real parts of error function with s = 0 (the 

left) and s = 1 (the right) for y 2 with uj = 1000. The fourth row: the real 

parts of error function with s = 2 (the left) and s = 3 (the right)for uo with 
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Figure 4.8: The top row: the real parts of error function with s = 0 (the 

left) and s = 1 (the right) for with oj = 100. The middle row: the real 

parts of error function with s = 2 (the left) and s = 3 (the right) for with 
Lo = 100. The third row: the real parts of error function with s = 0 (the 

left) and s = 1 (the right) for y^ with lo = 1000. The fourth row: the real 

parts of error function with s = 2 (the left) and s = 3 (the right)for with 
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Figure 4.9: The top row: the real parts of error function with s = 0 (the 

left) and s = 1 (the right) for 1/4 with oj = 100. The middle row: the real 

parts of error function with s = 2 (the left) and s = 3 (the right) for 1/4 with 
Lo = 100. The third row: the real parts of error function with s = 0 (the 

left) and s = 1 (the right) for ^4 with lo = 1000. The fourth row: the real 

parts of error function with s = 2 (the left) and s = 3 (the right)for va with 
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Figure 4.10: The top row: the real parts of error function with s = 0 (the 
left) and s = 1 (the right) for with oj = 100. The middle row: the real 
parts of error function with s = 2 (the left) and s = 3 (the right) for with 
Lo = 100. The third row: the real parts of error function with s = 0 (the 
left) and s = 1 (the right) for y^ with lo = 1000. The fourth row: the real 
parts of error function with s = 2 (the left) and s = 3 (the right)for 1/5 with 
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